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Abstract 

The lowest energy configurations for N equal charged particles 
confined to a thin conducting disc have been investigated in detail up 
to = 160 and in outline for further values up to = 500. For 
all values of N up to 160 the particle configurations can be described 
in terms of concentric shells. The number of perimeter particles p 
appears to be simply related to and to the mean radius of the out- 
ermost internal shell. Justification for these relations is obtained from 
a simple model based on the well-known distribution of continuous 
charge on a conducting disc. 
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1 Introduction 



If N identical charged particles are confined to the interior of a conducting 
hollow sphere, all the particles take up positions on the sphere itself - there 
is no stable, static equilibrium involving any particle with r < 1 [1]. The 
problem of how the particles arrange themselves on the sphere was originally 
posed by Thomson [2] and is still of great interest [see for example 3,4]. The 
corresponding problem for similar particles confined to a thin conducting 
disc has a more recent history. Given that in the sphere case all the particles 
must reside on the boundary, it is perhaps tempting to assume that in the 
2-D problem the N particles will take up equally-spaced positions around 
the perimeter. For small values of A^, this is exactly what is observed, but in 
1985 Berezin [5] reported the 'unexpected' result that for the lowest energy 
configurations with N > 11, the particles do not all reside on the perimeter 
- it is energetically favourable for some of the particles to be positioned in 
the interior. It was quickly pointed out, however, [6,7] that for very high 
N the configuration of discrete charges should approach the distribution for 
continuous charge, for which the result is well known and does indeed feature 
a non-zero charge density at all radii. 

For continuous charge, the charge density a{r) is calculated with great 
economy by Friedberg [8]. The equilibrium distribution of charge across 
the disc is shown to be equivalent to that found by projecting the uniform 
distribution on a sphere onto its equatorial plane. The result, also derived 
by conventional methods [9] is 

2a(0) 

'^(^^ = (1 -,2)1/2 (1) 

where a(0) is the central charge density. The charge density becomes infinite 
at r = 1. 

For discrete charges, at least for low values of N, the particles are arranged 
in a system of concentric shells. Detailed calculations for N < SO have been 
performed by Nurmela [10], who reviews earlier work and has reported the 
energies for each value of at what is thought to be its optimal configuration. 
Oymak and Erkoc [1 1] have also reported configurations and energies for < 
90. There is also a body of related work [see for example 12,13] in which the 
particles are confined to the disc by a parabolic potential {V{r) = l/2mwQ^r'^) 
rather than by hard-wall confinement (r < l,V{r) = 0; r > l,V{r) — oo). 
Here we use hard- wall confinement only and study the range 21 < A^ < 160 
in detail and 170 < N < 500 in outline, obtaining a relationship between the 
number p of perimeter particles and N. 



2 Methods 

2.1 Basic Methods 

All the programs used in this work employ the same relaxation routine. An 
initial distribution of particles is generated within the unit circle and the 
distances between all pairs calculated. The net Coulomb force on the ith 
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particle is 



(2) 



and each particle is then moved, subject to remaining within the perimeter, 
in the direction of the net force and by a distance chosen to be equal to 
F/N"^. The procedure is then repeated and periodically the total energy 



is calculated. In the simplest version of the program, execution is halted 
when the total energy fails to improve by a fixed amount (for example, 
AE = 0.0001) over 1000 iterations. As N increases the calculation of the 
intcrparticle forces eventually dominates the processing, requiring a time 
0(A^^), so in some versions of the program there are optional survey modes 
in which AE is increased to 0.001 or even 0.01. 

The system of N particles may encounter a local energy minimum from 
which no further progress may be made. Fig.l shows part of the spectrum 
of such metastable states for = 100, organised by the number of particles 
on the perimeter p{N). These alternative endpoints can account for a large 
fraction of the computing resources if no steps are taken to limit them. The 
approach taken here has been first to restrict p{N) to a specified range, but 
only after an initial survey has established the most likely p value. Some 
versions of the program also employ a simple form of annealing in which the 
positions of the interior particles are repeatedly disrupted with the aim of 
avoiding local energy minima. 

For programs running with a specified value of p, the required number of 
particles are set uniformly around the perimeter and the remaining N — p 
at random in the interior. If no particular value of p is specified, however, 
care must be taken with the initial distribution because the perimeter count 
tends to freeze early in the calculation. As Nurmela [10] points out, particles 
arriving first at the perimeter tend to inhibit further particles from joining 
them. If the initial distribution is totally random, the final states tend to 
be drawn from those on the left of spectra such as Fig. 1. If all particles 
are placed initially on the perimeter, the final states are biased to those on 
the right of the diagram. To achieve a range of p values centred close to 
the optimal value, we have employed initial radial distributions of the form 
r — 1 — A * random where the factor A ~ 0.1 must be adjusted slightly for A^. 

2.2 Three Stages of Computation 

The initial stage of computation was concerned only in the preparation of a 
list of candidate values of p{N). For 21 < A^ < 80, these were taken from the 
well-established results of [10]. For 81 < A^ < 160, various early versions of 
the program described below were used. Support for the full list of candidates 
was then obtained by generating an energy spectrum of final states similar 
to Fig.l and for each A^. Energy determinations of the highest accuracy were 
not sought and a basic program without preselection of p and terminating at 
AE = 0.01 was used. In every case the envelope of the lowest energies was 



N ^ 
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Fig.l. Part of the spectrum of final states for N = 100, organised 
by value. The parabola traces the envelope of the lowest energies 
for each p and defines the zero on the energy scale. 



parabolic in p and was fitted to E = -Emin + a{p — Pmm) ■ Here, -Emin is a 
lower energy than the optimal value associated with an integer p-value, and 
plays no further part in this work. The values of pj^m are included in Table 
1 and the variation of a with is shown in Fig.2. 

The candidate integer values of p were then used in the main stage of com- 
putation. For each in the range 21 < iV < 160 a calculation was performed 
with a predetermined perimeter count set in turn to p — 1, p and p + 1. For 
each p-value, an initial distribution was generated which was then allowed 
to evolve during 100 sessions each of 10000 iterations. Between sessions, the 
position of each internal particle was modified by a random vector with am- 
plitude s * random. At each annealing stage, the value of s was reduced by 
1% of its initial value, which was taken as 0.4. This choice meant that in 
the initial stages for most A^ values, several of the internal particles moved 
outside the perimeter, forcing the generation of a new distribution. This 
annealing procedure is clearly less sophisticated than the formal simulated 
annealing algorithm first applied to this problem by Wille and Vennik [14], 
which requires that each particle is moved in turn and the new configuration 
is then accepted or not on the basis of the change in energy. 

The confirmation stage was a repeat of the main stage. The only data 
carried forward from the main stage were the p- values of the best results, but 
this time the variants with p — 1 and p + 1 were not calculated. The confir- 
mation stage was performed twice, and the two sets of results were compared 
with the results of the previous stage. The lowest energy of the three results 
for each A^ is shown in Table 1. A further run was then performed with an 
alternative form of annealing - instead of s = 0.4, s was taken as 1 — r, where 
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Fig. 2. Variation of the parabolic fit parameter a with N for 
21 < < 160 and for selected values up to = 500. The typical 
standard deviation for repeated measurements on a value of a is 
0.01. 



r was the radial distance of each internal particle from the centre. This pro- 
duces a less extreme disturbance of the particle positions, but was not found 
to offer any advantage. Though the bulk of the results were identical with 
those already obtained, 21 had marginally higher energies and only for two 
values of was a small improvement (one digit in the final quoted decimal 
place) obtained. 

2.3 Results for 21 < iV < 160 

For all 21 < < 80 the lowest energy was produced at the central p- value 
tested, and the energies and configurations agree in every particular with 
those of [10]. The results in this range are listed here for completeness and 
to include the measurements of pmin and of the mean radius of the outer- 
most inner shell r2. In the range 81 < < 160, for all but one value of 
the lowest energy found in the main stage was again produced by the cen- 
tral p-value of the three tested, and the energy determined was equal to or 
slightly below (one or two places in the final figure) the value noted at the 
survey stage. The small reductions reflect the superior performance of the 
100 restart/annealing stages compared with earlier versions of the program. 

2.4 Results for N > 160 

For N=170,180...290 the spectrum program was run without further investi- 
gation, and the values of a and Pmin only are shown in Table 1. The parabola 
constant a (Fig. 2) is well-behaved, and appears to settle to a value of around 
0.36 as A^ is increased, suggesting that the envelopes of all spectra eventually 
have a similar shape. This predictability is supported by the good agree- 
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ment between Pmin and the optimal value ofp{N) for N < 160 (Table 1) and 
suggests that only a limited range of p-values near p^nm need be investigated 
if computing resources are at a premium. Accordingly, for = 200, 300, 
400 and 500 the spectrum program was first used to determine Pmin and the 
annealing program was then run for one or two candidate values of p. Full 
configurations were not obtained for these values and the listed energies are 
not necessarily optimal. 

3 Shell Structure 

Three points were noted: 

(1) In the range N < 80, the configuration for + 1 is, with one exception, 
the same as that for A" with the addition of one particle to an existing shell 

or the creation of a new shell with a particle placed at the centre. The single 
exception is the pair A^ = 55 (5-13-37) and A^ = 56 (1-6-12-37). For the 
range 81 < A^ < 160 there are four exceptions to the usual pattern: A^ = 
97/98, 117/118, 150/151 and 152/153. 

(2) Oymak and Erkog [11] have suggested that new shells are created when 
A^ is given by 

m 

A^ = 12 + 5^(10 + 8n) (4) 

n=l 

This is equivalent to N=(2m+3)(2m+4) and suggests that new shells are 
created at A^ = 12, 30, 56, 90, 132, 182.... The results obtained here, however, 
show the next new shell after A^ = 90 appearing at A^ = 134. 

(3) Of the five possible 2D Bravais lattices, it is the triangular lattice that 
produces the lowest energy packing for identical charged particles [15], and 
this arrangement is observed in the central region for many values of A^. 
Nevertheless, for A^ < 90 previous authors have also described the optimal 
configurations in terms of concentric shells, and this property appears to 
be maintained up to A^ = 160. Beyond this value, the first evidence of a 
departure from this behaviour may occur at A^ = 185 (Fig. 3). Repeated 
attempts have been made to improve upon this structure and its energy 
{E = 23129.9182) but no lower energy state has been found. A^ = 185 may 
be an initial, isolated example of a later trend, because preliminary results for 
the remaining A" in the range 161 < N < 200 all appear to have well-ordered 
shell structures. 
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Table 1. Results for 21 < N < 72 
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Table 1 (cont.) Results for 73 < N < 124 
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Table 1 (cont.) Results for 125 < N < 160 and selected values 

N < 500 
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Fig. 3. The lowest energy configuration for N = 185, the first A^- 
value for which a simple concentric shell description appears to be 
inappropriate. The configuration has p = 87. 



4 Total Energy 

Erkoc and Oymak [16] have investigated the total energy, using the expression 

E{N) = ^N^ + h2N^''^ + h^^N (5) 

The leading term is the exact expression for E{N) as N —>■ oo. The remain- 
ing terms are the first two of an infinite series representing the correction 
for finite N. Erkog and Oymak fitted their data in the range 2 < N < 90, 
obtaining 62 = —1.5599728 and 63 = 0.9509338. Repeating the analysis 
with 2 < < 160, we obtain 62 = -1.5593651 and 63 = 0.94290006. These 
values predict E{500) = 179386.7591 compared with the observed upper 
limit of £"(500) = 179375.07. If a fit to a limited range of the data with 
100 < < 160 is performed, we obtain 62 = -1.5611781 and 63 = 0.9614771. 
These values predict _E(500) = 179375.7771, significantly closer to our best 
result for = 500. 

5 Perimeter Particles 

If the shell structure is indeed lost as A^ increases, it nevertheless seems likely 
that the perimeter particle count p{N) will survive as a useful parameter 
and should be recognised as an integer sequence [17]. As support for such 
a sequence, an attempt has been made to express p as a function of A^ up 
to and beyond A^ = 160. This has been achieved as a by-product of an 
investigation into a faster method of identifying optimal configurations for 
higher values of A^. 

The increase in processing time required to obtain reliable solutions at 
higher A^ suggests that some pre-preparation of the particle distributions 
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may be advantageous. Artificial distributions on their own are unlikely to 
deliver optimal solutions which are not, for the most part, made up of perfect 
circles and in any case the shell structure is unlikely to survive beyond some 
value of N. Nevertheless, it was felt that such an approach could extend the 
range of reported results if suitable patterns of shells were used as inputs to 
the procedures described above. If a sufficient number of different patterns 
are used as seeds, such that at least one is not separated from the optimal 
final state by a potential barrier, the need for the repeated annealing stages 
may be avoided. 

Table 1 and [10,11] all show 4-10-18-48 as the optimal configuration for 

N — 80. To generate such a sequence together with a number of nearby al- 
ternatives, the approach adopted here is to begin with charges on the disc 
in an equilibrium distribution of continuous charge. We then regard the tran- 
sition from such a distribution to a lowest-energy configuration of N discrete 
charges as a process of condensation. To arrive at the correct final configu- 
ration it is necessary to establish the boundaries separating each unit charge 
on the continuous distribution. (An alternative method of establishing such 
boundaries, the Voronoi construction [18] has been applied to this problem 
by Bedanov and Peeters [13] who report results for hard-wall confinement 
with A^ = 50 and A^ = 230 only, obtaining a configuration of 1-5-13-31 for 
A^ = 50.) After integration of Eqn.l, the fraction of charge outside radius 
r is given as /(r) = (1 — r'^y/'^ and between any two radii and rb with 
'"b > "/"a the fraction of the charge is g = (1 — r^'^y^'^ — (1 — rh^Y^"^. If such 
an area is divided into n equal sectors such that each contains exactly the 
amount of charge to represent one particle, then g/n— 1/N. We begin at the 
perimeter, with all A^ particles still available for allocation and = 1, and 
ask how many particles n should be assigned to the outermost shell. Clearly 
at this stage the value of n may fall in the range 1 < n < A^. If n = A^ is 
chosen then = 0, but for all other values we obtain < ra < 1 (Fig.4). 
For all possible values of n we calculate the length of the boundary of the 
cell representing one particle and accept the value of n giving the smallest 
boundary. This defines p = n. The corresponding value of ra is then taken as 
the outer radius contributing charge to the next shell, the remaining number 
of particles is adjusted to A^ — n and the process is repeated until no particles 
remain. 

Applying the above procedure for A^ = 80, we obtain 2-8-14-22-34. For 
A^ = 160 the configuration (Table 1) is 3-9-15-22-33-78, and the procedure 
yields 1-7-13-20-28-37-54. But it was noted that by making a simple modi- 
fication to the calculation of the boundary, initially for the first (perimeter) 
shell only, much closer agreement could be obtained. A weighting factor 
w was introduced for the radial part of the boundary, so that its effective 
circumference s was taken as 

s = 2w{rb - Ta) + 27r(rb + r^)/n (6) 

and the minimum value of s is used to select ra. With wi — 0.3 for the 
first shell and W2,W3... — 1.0 for all subsequent shells, the new configuration 
obtained for A^ = 80 is 4-10-17-49, and for A^ = 160, 3-9-15-23-31-79. 

A fuller investigation was then undertaken to establish the values of the 
weights that generate the exact optimal configuration for each A"- value in the 
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Fig. 4. Working inwards from the perimeter, various numbers of 
particles (initially 1 < n < A^) are allocated to each shell such that 
for each choice of n, each cell contains a single charge, using the 
distribution given by Eqn.l. The number adopted is the one pro- 
ducing the lowest cell perimeter (shown in blue), taking into ac- 
count the weighting factor which multiplies the radial component. 
The procedure is then repeated for the next inner shell until no 
particles remain. 



range 20 < < 160. At a particular value of A^, wi for example can take a 
range of values and still deliver the correct p, so each produces a maximum 
and a minimum value for each weight. The results for the perimeter weight 
Wi are shown in Fig. 5, which also includes the results for selected values with 
A^ > 160 obtained from A^ and p only. It is found that as A^ increases, the 
correct perimeter count p may be obtained with a value of Wi close to 0.30. To 
generate a collection of candidate configurations for further processing by the 
relaxation method, it is only necessary to repeat the condensation routine 
with weights wi,W2--- each varying in a narrow range about its adopted 
central value. Further results beyond A^ = 160 using this method will be 
presented in future work. 

The eventual similarity of the results for wi suggests that the consequences 
of a fixed value of this weight should be investigated. We assume that there 
is some radius r such that the fraction of charge outside r accounts for the 
perimeter particles, and write p = N{1 — r^)^/^ = A^(l — r)^^^(l + r)"^''^. 
Substituting this version of p into Eqn.6 and setting n = p, ra, = r and 
Tb = 1, we obtain 




Hence, to select the minimum boundary length 




1/2 



-2wi + 
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0.45 




0.25 

0.2 I ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

100 200 300 400 500 600 
Total Particle Count N 

Fig. 5. The range of values of wi that generate the correct con- 
figuration for 20 < < 160 and selected values up to = 500. 
The results for > 160 are obtained without knowledge of the 
full configuration and use and p to measure wi only. 



= for p(l -r) = — (7) 

Wi 

If we take Wi = 0.3 and identify r as the radius of the outermost inner 
shell r2, minimising the boundary leads to the relation p{l — = 10.47. 
Fig. 6 shows the variation of p{l — with for all the results shown 
in Table 1. The mean value in the range 100 < N < 160 is 11.461, and 
the higher values of remain close to this value. The quantity p(l — T2) 
appears to tend to a constant value, though we note that r refers to the inner 
boundary of the region producing the perimeter particles, and hence the outer 
boundary contributing the second shell. This radius is not necessarily the 
radial position of the second shell particles r2, so the values of the constants 
need not agree. 

Finally, substituting for the remaining r in Eqn. 7 produces 

For large N such that p/N <^ 1 the square root term may be replaced by 
its three leading terms, resulting in a quintic equation in p. The two leading 
terms of the solution are 




If we again take Wi = 0.3, Eqn. 8 becomes p = 2.76N^^^ — 1.75. Using the 
results for all N listed in Table 1, Fig. 7 shows the dependence of pmin upon 
A^^/^, yielding a best fit of p = 2.7922iV2/3 - 3.723 ± 0.021. Once again. 
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Fig. 6. The relationship between p and r2 for 21 < < 160 and 
for selected values of up to = 500. 



the results are not strictly comparable - both use the same p and N data, 
but Eqn.8 relies upon our assumption that p/N <C 1. Nevertheless, the 
simple method of charge partitioning using only the single parameter wi as 
an input produces a relationship between p and which appears realistic. 
The estimated error on the fitted constant term was then used with the 
coordinates of the centroid of our data (Ar2/3 = 21.7988, p^^^ = 57.1445) to 
predict values of Pmin beyond A^ = 500. We obtain pmin = (3.723±0.021)s + t 
where s = 0.04587A^^/^ — 1 and t = 2.6215A^^/^. Although such extrapolations 
must clearly be treated with caution, we have used this to estimate the 
range of Pmm for A^ = 1000, obtaining pmin = 275.50 ± 0.08. Efforts with 
the spectrum program have so far yielded a best energy oi E = 736985.16, 
obtained at p = 276. This p-value must be regarded as preliminary and the 
energy as an upper limit, but the relation between p and A^ does appear 
to be quite robust. If p(lOOO) = 276, this requires (Fig.5) Wi in the range 
0.292-0.294. 

Atiyah and Sutcliffe [19] have designated the A^-values at which the latest 
particle is added to the interior of the circle rather than to the perimeter as 
'jumping points'. The sequence of such points begins A^ = 12,17,19,22... 
and in our notation these are the points of unchanging p. This suggested 
an alternative method of fitting Eqn. 8 by searching for ci and C2 such that 
the nearest integer to ciN'^^^ — C2 gave the correct value of p in the highest 
possible number of cases. It was found that in the range 21 < A^ < 160 all 
but eight values of p could be correctly generated with ci = 2.7777 and 
C2 = 3.3203. If a 100% success rate had been achieved, the jumping points 
would have been directly calculable from the expression for the p-values. 
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Fig. 7. The variation of Pmin with N"^^^. The typical standard 
deviation on a measurement of Pmin is 0.04. 



6 Conclusions 

We have obtained solutions for the optimal distributions of A'^ charged parti- 
cles on a unit disc, extending the known configurations to = 160. Although 
it appears that well-behaved patterns of concentric shells may become less 
prevalent after A^ ~ 185, it is suggested that the perimeter count p{N) should 
survive as a useful parameter, and tentative values for pmin have been ob- 
tained for further selected values of A^ up to 500. Over the range studied 
p varies according to p = 2.7922A^2/3 _ 3 723. if this expression is used to 
determine p(lOOO) the result appears to be supported by a preliminary calcu- 
lation. A simple procedure for assigning particles to concentric shells appears 
to support the observed relationships between p and A^ and between p and 
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